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We prove the existence of chaos in temperature in the Sherringhton-Kirkpatrick model. The effect 
is exceedingly small, namely of the ninth order in perturbation theory. The equations describing two 
systems at different temperatures constrained to have a fixed overlap are studied analytically and 
numerically, yielding information about the behaviour of the overlap distribution function Pt-^^t2{<i) 
in finite-size systems. 
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I. INTRODUCTION 

In this work we address the question of chaos in temperature in the Sherringhton-Kirkpatrick (SK) model. According 
to the Parisi solution Q the SK model has many equilibrium states at a given temperature. These states are correlated: 
the overlap between different states can take a whole range of values from zero to the self-overlap qea- A natural and 
simple question is: what is the overlap between states at different temperatures? According to the chaos in temperature 
hypothesis the overlap between any couples of equilibrium states at different temperatures is zero. Consequently, the 
probability distribution of the overlap of configurations equilibrated at different temperatures P{q) should be a delta 
function centered at g = 0. One can be more quantitative considering the convergence to a delta function of the 
function PAr(g) corresponding to a finite system size N. Standard arguments tell us that for large N this quantity is 
given by the large-deviations formula 

PN{q)^eM-NAF{q)] (1) 

the quantity AF{q) — F{q,Ti,T2) — {F(Ti) -\- F{T2)) is the free-energy cost, i.e. the difference between the free 
energy of a system of two copies at different temperatures of a spin-glass with the same realization of the disorder 
constrained to have a given value q of the overlap, and the sum of the free energies at temperatures Ti and T2 (actually 
here and in the following we call "free energy" the opposite of the logarithm of the partition function, without the 
usual /3 factors). Consequently we can talk of a small or large chaos effect depending on the function AF(g, Ti,T2); if 
AF{q, Ti, T2) is small the convergence with TV of the Pat (9) will be slower and the chaos effect will be seen only when 
considering large-size systems. In particular if AF{q) turns out to be zero for a finite range of values of the constraint 
there is no extensive chaos; furthermore a null AF is expected if the P{q) is definitely different from a delta centered 
on zero. In the following we shall compute the quantity AF(q). 

In a recent paper , one of us reexamined the problem through a perturbative expansion in the replica framework, 
showing that chaos, if at all present, should be a very small effect. Actually, this is a quite old problem H, ^, ||J^, |j, 
and has gained new attention in connection with the experimental effects of memory and rejuvenation (see ToTand 
references therein). Recently some interesting results have been found for other models. Sales and Yoshino pT| ] 
have convincingly shown that there is chaos in the directed polymer in random media (DPRM) in 1+1 dimensions 
due to the interplay between energy and entropy, a mechanism which had been postulated in the droplet theory for 
finite dimensional spin-glass ||ll| . On the other hand the existence of high correlations between systems at different 
temperatures has been proved in a well known class of mean-field spin-glass models [|l^ . Some evidence against chaos 
in the naive TAP equations has been found in ||l3|, a result which has been questioned in iQ. Powerful numerical 
efforts have been devoted by Billoire and Marinari (BM) to investigate chaos in temperature in the SK model ||l^, ^ ; 
even though they don't actually see the effect, there is some evidence that it may be extremely weak, thus observable 
only on very large lattices. 
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II. THE VARIATIONAL EQUATIONS 

In the replica framework one obtains a variational expression for the free-energy functional F{qc,Ti,T2) of two 
systems with a constraint Qc ||, 0, We have expressed F(q(., 7i, T2) in terms of a Parisi differential equation. The 
variational parameters are three n x n matrices Qi = (0, gi(a;)), Q2 = {0,q2{x)), and P = {pj^,p{x)): 

F{qc,T^,T2) = ^^^^ + ^(- [\f(:c)dx + 2q,{l))+^(- (\l{x)dx + 2q2{l 



4 4 4 V 7o 7 4 

^\\x)dx\ - - ,.)^ - /(0,0,0) (2) 



The function /(a;, yi, 1/2) is defined through the following differential equation: 

df qi fd^f , fdfY\ , 92 fd^f , (Of 



dx 2 \ dy\ ' \ dyi J j 2 I Qy^ V ^2/2 
\0y10y2 oyidy2j 

With initial condition: 

/(1, 2/1,2/2) = ln[2e^ cosh[/3i2/i + /322/2] + 2e-* cosh[/3i2/i - /?22/2]] (4) 

<5 = /3i/32(pd-p(l)) (5) 

Notice that the initial condition depends on the difference pd ^p{l). In order to maximize the free energy functional 
Lagrange multipliers P(x,2/i,2/2) are used much as in the standard single system problem [ p^ , pj| , then introducing 
the following quantities 

mi[x,yi,y2) = —-^ , ^2 (a;, 2/1, 2/2 = 73-;^ (6 
Pi 02/1 P2 oy2 



and the operator 



,dyidy2 Ti dy2 T2 dyi 

the saddle-point (SP) equations for the constrained systems read: 

dP qi fd^P 2x d \ q2 fd^P 2x3 

- (Pmi ) + — — 2- - — — (Pto2 ) 



d qi f 2x d \ q2 f d"^ 2x d 

. f d'^ X d X d \ 
+ P o — o ^ ^^"1 T 1^ ^"''2t; — , (7) 



aa; 2 \dyl T, dyi"^ ''J 2 \dyl T2 dy2 
( P X X \ 

\oy1OV2 Tidy2 T2 dyi ) 

^mx = ; ^■m2 = (9) 

qi{x) = J dyidy2Pmj; q2{x) = J dyidy2Pml; p{x) ^ J dyidy2Pmim2; (10) 

The initial conditions for the functions mi (x, 2/1, 2/2), "^2(2;, 2/1, 2/2) can be obtained from and (||); for P(a;, 2/1, 2/2) 
we find 

P(0,2/i,2/2) = <5(2/i)^(2/2) . (11) 
The presence of a constraint q^ leads to the following equation 

f, , n.-, ^e'^cosh[/3l2/l -h/32y2] -e-'^cosh[/3i2/i -/322/2] 

9c = / dyidy2P(l,yi,y2)— tj^ —5 — r- — ^5 rj^ » — T • (^2) 

7 e"' cosh /3i2/i + /322/2 + e * cosh /3i2/i - /322/2 



If (5 = the initial condition (|J) separates and eq. ( |12[ ) implies p(l) = Qc- Finally, we recall that the derivative of the 
free energy with respect to the constraint can be expressed as 

OF 

e = — ^ PiP2{Pd ~ Qc) ■ (13) 
oqc 
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III. THE ^1 = /?2 CASE 

When Ti = T2 there are two possible situations depending on the value of the constraint Qc- For all values of qc 
inside the support of the overlap distribution function P{q) we have solutions with 6 = and therefore e = 0. As a 
consequence for < Qc < Qea we have /S.F = as it is to be expected since the function P{q) is non trivial. The 
solutions of the SP equations satisfy pa = 9c and are: 

qi{x) = q2{x) = p{x) = qparisii'^.x) < X < \xparisi{j)d) 

qi{x) = q2{x) = p{x) = Pd \xparisi{Pd) <X < Xparisi{Pd) 
qi{x) = q2{x) = qparisi{x);p{x) = Pd XparisiiPd) <X<1 (14) 

These solutions were first proposed in for the truncated model. In it has been proven that they exist in 
any models with Replica-symmetry breaking (RSB) by noticing that they are a permutation of the standard Parisi 
solution. 

When lies outside the support of the function P{q) there are off-equilibrium solutions of the SP equations 
with (5 7^ and qi{x) = q2{x) = p{x) These solutions can be obtained by solving the SP equations (||j9|jl^) 

perturbatively in Sq ^ qc — qEA{T) and r = T — Tc- This has been done at the lower orders in where it was 
found that near the critical temperature the free-energy difference grows with 6q'^ in the off-equilibrium region. We 
pushed the computation to very high orders and we can safely claim that there are no 6q^ terms in the expansion. 



Recently this result was confirmed numerically |20 . The off-equilibrium phase can be described also by considering 
what happens for — 1 rather then for qc ~ qsA- In the following section we will discuss this limit which can also 
be used to obtain another series expansion. 

IV. THE qc = 1 CASE 

In this case the only possible couples of configurations are those in which a configuration is coupled with itself. 
It is easy to realize that the thermodynamics of such a constrained system is equal to that of a single system at 
temperature T12 = l/(/3i + ^2)- It is also natural that qi{x) = 92 (a;) = p{x). By looking at equation (|lj) we see that 

= 1 requires 6 = 00. In this case the initial conditions on mi and m2 read 

mi(l, 2/1,1/2) = m2{l,yi,y2) = tanh[/3i2/i 4-/322/2] (15) 
By making the change of variables 

S = r, .n — ' D^y2-yi , (16) 

Pi + P2 

the SP equations (^,^,Q reduce to those of a single system at temperature T12, as expected. The case Ti = T2 
has been recently discussed in It is interesting to notice that the constrained model with Qc = 1 has a critical 
temperature whose value is exactly twice as the critical temperature of the free model. Therefore if we consider a couple 
of constrained systems such that Ti > Tc, T2 > Tc but T12 < Tc we shall have a phase transition from a paramagnetic 
replica-symmetric (RS) phase to a RSB phase by changing the value of the constraint qc- For completeness we report 
the paramagnetic RS solution of the SP equations: 

qi{x) = q2{x) = p{x) =0 pd = -^-T- arctanh qc (17) 

P1P2 

Notice that pd is divergent for gc — 1; therefore, according to (|l^), the derivative of the free-energy is divergent too. 
This result is valid also in the RSB phase. For Ti = T2 = T, 1 < T < 2 the critical value of the constraint qc turns 
out to be gcrit(T) = T — 1. For qc near (jcrit the order parameter is small and we solved the equations perturbatively 
in qc — gcrit- We also obtained a perturbative expansion for Ti ~ 2, T2 ~ 2 and qc ~ 1. Indeed, in this region the 
three functions q\[x), q2ix) and p(x) are small since the value (/c = 1 corresponds to a free system at a temperature 
near Tc = 1. 

V. THE /3i / /32 CASE 

In 1^ an analytical picture for the absence of chaos in spin-glasses was proposed based on the possible existence of 
specific solutions of the SP equations (0,0,^01). The structure of these solutions is depicted in figure (IT]): in the small- a; 
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FIG. 1: The solutions discussed in The real solutions are expected to pick up small corrections induced by chaos, 
actual difference between the three functions in the small- a; region is much smaller than in the figure. 



The 



region the three functions qi{x), q2{x) and p{x) are all different till they reach the point Xc where p{x) = Pd = Qc] 
then for x greater than Xc p{x) remains constant while qi{x) and 92(2;) after an intermediate plateau are connected 
continuously to the corresponding free Parisi solutions. These solutions may exist for values of the constraint from 
zero to a maximum value where the two plateaus of the function at the higher temperature merge (Ti > T2 in figure 
(|l|)). These solutions reduce to eq. in the limit Ti — > T2 and are intrinsically non-chaotic due to the relation 

5 — [||. However, their existence must be checked explicitly by attempting to solve the SP equations (^,||,^0|) with 
the ansatz 5 — 0. In |^ this has been checked up to the fifth order in a perturbative expansion near the critical 
temperature, i.e. it has been found that Ai^ is a quantity smaller than the fifth order. 

Recently it has been proven that these solutions exist in the spherical spin-glass models with multi-p spin interaction 
[ p2[ . In that context it has been shown that they imply ultrametricity between the equilibrium states at different 
temperatures so that the Parisi tree is essentially the same at all temperatures. Furthermore, arguments have been 
advanced in order to show that the reference free energies of the clusters of states do not change with temperature, 
much as in the generalized random-energy model. Perfect ultrametricity is a special property of those models due to 
the fact that the functions qi{x), q2{x) and p{x) are all equal in the small- a; region. In general, however, we believe 
that the physical interpretation of the solutions proposed in is similar. Indeed, we think that the non-chaotic 
solutions proposed in yj are the mathematical implementation of the old idea of a unique Parisi tree which bifurcates 
when lowering the temperature plus the additional hypothesis that the ordering of the reference free energies of the 
clusters of states does not change when lowering the temperature. It would be interesting to check whether the logical 
equivalence between these solutions and the previous picture can be established explicitly through the cavity method 
as was done for the Parisi solution fl| . 



A. Perturbative Expansion 

In order to deal with the SP equations of the SK model we resorted to a high-order computer-assisted 

perturbative expansion using the methods we recently applied to the Parisi solution fl^ . The procedure is conceptually 
straightforward: one tries to solve the SP equations with the non-chaotic ansatz of figure (|l|) assuming that near the 
critical temperature the functions (j'i(x), 92(2;), p{x) are small and the x regions where they vary are also small. 
Therefore one expands all the quantities of interest in series of y, a;, ti, T2 and qc, {Ti — 1 — Ti); the first is of order 
1/2, the others are of order 1; the function P{x, y), which is singular, can be treated considering its y cumulants. This 
is essentially the same computation that was carried on in Q by hand up to the fifth order of the free-energy. It was 
found that up to that relatively high order the SP equations admit of a non-chaotic solution with (5 = for values of 
qc between and some maximum value near the self-overlap of the states at the higher temperature. Instead we have 
found that at ninth order in the free energy it is impossible to find a solution of the SP equation with 5 = 0: 

the non-chaotic solutions don't exist in the SK model. Accordingly there must be a free-energy cost in imposing the 
constraint, which one expects to be of ninth order. 

Since off-equilibrium solutions must have 5^0, things become very complicated; indeed, as explained in [||, when 
5 the three functions qi{x), q2{x) and p{x) are coupled in the whole small- a: region, i.e. if one of them is varying 
in some interval of x the others must vary too. We guess that off-equilibrium solutions have the same structure of 
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the non chaotic solutions of figure (|l|) plus small corrections of order S, which is expected to be of the seventh order. 
Therefore we have tried to solve the SP equations (^,| 10) with a complicated ansatz made up of five x regions on 



each of whom the functions qi{x), q2{x) and p{x) have a different expansion. Unfortunatel y, i n spite of the very high 
number of parameters of this ansatz, we have not been able to solve the SP equations (^|Jg|Jl0|) explicitly . Therefore, 
in order to estimate the free-energy cost we plugged the same ansatz directly into equations (|2| ^) and extremized the 
free energy with respect to the variational parameters. We obtained the following estimate for the free-energy cost of 
constraining two copies of a SK spin-glass at temperature Ti = 1 — ti and T2 = 1 — T2 to have a fixed overlap qc- 

AF^I |g,r(ri-r2)2 (18) 

As expected, this quantity turns out to be of the ninth order. Let us comment that although in principle this must be 
considered an upper bound for the real value, it is most likely exact. Indeed, in we observed that the order of the 
function q{x) relevant for the exact value of the free energy is much lower than the one expected by power counting 
arguments; for instance if we maximize the free energy at order n, we are able to determine the order parameter at 
an order near n/2 rather than the order n — 2 obtained by power counting. Therefore even if we know that our ansatz 
is unable to reproduce the correct solution at the seventh order, most likely this has no influence on the value of the 
free energy at the ninth order. 




FIG. 2: The functions qi{x), q2{x) and p(x) obtained numerically for qc = 0.3, Ti — 0.5 and T2 = 0.4. The dotted lines 
represent the corresponding free solutions. Within the numerical precision they coincide with those of figure (|l|); the accuracy 
could be raised by the methods discussed in 



B. Numerical Integration 



We also integrated the SP equations numerically via a pseudo-spectral code as we did in Q for the standard Parisi 
solution. Since the fields have three variables, time and memory required by the computation grow and the obtainable 
precision decreases, preventing us from seeing the chaos effect. The advantage of this method, however, is that we 
can investigate the form of the solutions in a non-perturbative way. To obtain the three functions qi{x), q2{x) and 
p{x), we applied the iterative scheme of |^ to the SP equations flpjiol) ; we imposed at each iteration p{x) < qc 
and we set by hand S — 0; indeed S is expected to be much smaller than the attainable precision. Consistently with 
this assumption, the free-energy cost computed integrating (|^) and (^) turned out to be zero within the precision. 
Remarkably the three functions qi{x), q2{x) and p{x) have precisely the structure proposed in |Q (cfr. figures (|l|) and 
(H)). The chaos effect determines small corrections which should be seen at higher precision. When qc is greater than 
some Qmax the plateau of the function p{x) stops at g,nax irrespective of Qc, and in order to obtain higher values of 
p{x) we should impose a finite value of 5. In the isothermal case this happens when one goes out of equilibrium for 
qc > qEA\ when Ti 7^ T2 this effect marks a change in the behaviour of Ai^ which is much larger for qc > (/max (see 
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discussion in the conclusions). We believe that the chaos effect is beyond the attainable numerical precision of our 
code at the present, while improvements for the functions qi{x), q2{x) and p{x) are at hand. 



VI. CONCLUSIONS 



We settled the longstanding problem of chaos in temperature in the SK model: chaos in temperature is present 
but it is an exceedingly small effect: AF is a ninth order quantity. The analytical picture for absence of chaos in 
temperature proposed in |^ doesn't apply to the SK model, while it holds in the multi-p spin spherical models [p^ ; 
this is an important difference between these two models whose behaviour has been often assumed to be similar. 

The effect is so weak that it can't be observed in present numerical simulations. Indeed, if we make an estimate 
for the function P/v(g) through the equations and ( p^ ) with parameters Ti = .6, T2 = .4 and N = 4096, which 



are the temperatures and size investigated in the latest numerical work by Billoire and Marinari (BM) |16 , we find 



a completely flat function; therefore we don't expect to see anything similar to the delta peak corresponding to the 
iV — > cxj limit. We notice that other corrections must be considered in order to explain those data; indeed, according 
to our results, a formula like -Pjv(9) oc exp[— A'^Ai^] in any case would give a -Pjv(<z) with an absolute maximum in 

Although we have not been able to determine exactly the off-equilibrium solutions corresponding to a non-zero 
constraint, we believe that they are similar to the non-chaotic solutions proposed in |^ plus small correction of order 
6 = 0{t'^); this is supported by the numerical integration of the SP equations (see fig. ||). As a consequence some 
of the features of the non-chaotic models discussed in are valid in the SK model at finite N . For instance, we 
expect that the support of the function Pniq) has a maximum Qmax located near the self-overlap q_EA(Tgrcat) of the 
system at the higher temperature. In the spherical models this happens because AF = for g < qmax and AF > 
for q > qmax] instead in the SK model (/max marks a crossover in the behaviour of the function AF which is ninth 
order for q < (/max while for q > gmax is much larger, most likely third order as in the isothermal case. This is clearly 
seen in the numerical simulations of BM and also confirms their hypothesis that the approach to a delta function 
doesn't proceed by a splitting of the peak at (/max but rather by the emergence of a third peak in q — 0. 
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